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EXECUTIVE  SUMMARY 


The  overall  goals  of  this  program  were  to  develop  fast  thermal  dissipation  mechanisms  for 
lightweight  polymer  matrix  fibrous  composites  based  on  moisture  evaporative  cooling, 
superconducting  heat  pipe  arrays,  and  thermal  ceramic  protective  coatings  for  future  Air  Force 
vehicles,  and  associated  propulsion  systems.  The  project  had  these  outcomes: 

•  Models  show  that 

o  Moisture  bearing  composites  can  delay  and  decrease  the  surface  temperature  rise  that  occurs 
when  exposed  to  90  °C;  however,  the  effect  is  limited  by  the  moisture  carried  by  the 
composite  and  the  time  taken  to  replenish  moisture  in  the  structure, 
o  Flowing  evaporative  cooling  has  a  similar  effect;  the  moisture  flows  through  the  structure  to 
cool  the  surface.  The  method  is  effective  for  low  relative  humidity  environments. 

•  Our  experiments  failed  to  produce  components  with  naked  internal  channels  for  fluid  circulation. 

This  outcome  suggests  that  channels  might  be  made  from  microtubes  that — when  combined  with  the 
reinforcing  fibers — remain  in  the  composite  to  direct  fluids  through  the  structure. 

•  The  superconducting  heat  pipes  that  the  proposal  was  based  upon  remain  controversial.  To  date, 
publications  are  few  and  show  inconsistent  behavior.  The  proposed  work  was  to  consider  the 
superconducting  heat  pipes  as  candidates  for  miniaturization — assuming  the  superconducting 
properties  were  functional  at  the  microscale.  Data  to  support  this  plan  is  unavailable  at  this  time. 
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OBJECTIVES 


•  The  overall  goals  of  this  program  are  to  develop  fast  thermal  dissipation  mechanisms  for 
lightweight  polymer  matrix  fibrous  composites  based  on  moisture  evaporative  cooling, 
superconducting  heat  pipe  arrays,  and  thermal  ceramic  protective  coatings  for  future  Air  Force 
vehicles,  and  associated  propulsion  systems.  The  thermal  ceramic  protective  coating  work 
stopped  in  early  2008  when  the  project  was  refocused  after  the  tragic  loss  of  the  originating  P.I., 
Dr.  Roger  J.  Morgan. 

•  Evaporative  moisture  rapid  surface  cooling  of  composites  will  be  evaluated  by:- 

o  Experimental  diffusion  of  moisture  from  high  moisture  bearing  fillers  embedded  in  a 
polymer  matrix,  in  conjunction  with  moisture  diffusion  modeling  studies 

o  A  network  of  microcapillaries  that  supply  water  to  the  surface  structure  from  an  internal 
reservoir. 

•  An  array  of  ultra-fast  conducting  heat  pipes  will  be  developed,  that  are  based  on  internal  thermal 
gaseous  increases,  which  rapidly  conduct  surface  heat  to  an  internal  heat  sink. 
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STATUS  OF  EFFORT 


The  project  stopped  technical  progress  in  December  2009  with  some  report  activity  continuing  into 
January  2010. 

MODELS 

Numerical  Evaluation  of  Active  Evaporation  Cooling  of  a  Composite  Panel  Using  the 
Finite  Element  Method 

Introduction 

Limiting  the  external  temperature  of  stealth  aircraft  has  become  a  major  concern  as  aircraft  have 
become  more  vulnerable  to  detection  due  to  recent  improvements  in  the  IR  signature  detection 
devices  of  missile  systems.  If  an  aircraft  deviates  from  its  surroundings  by  only  1°C,  detection  can 
be  achieved  at  military  useful  ranges  [1].  Current  stealth  aircraft  utilize  cooling  pipes  within  the 
wings  and  engine  exhaust  to  limit  external  temperatures  and  thereby  the  associated  IR  signatures. 
However,  more  efficient  and  instantaneous  cooling  mechanisms  are  needed. 

We  numerically  evaluate  the  potential  of  evaporative  cooling  of  water  as  an  alternative  mechanism 
for  the  cooling  of  stealth  aircraft.  Evaporative  cooling  of  liquid  water  occurs  when  the  surface  of  a 
body  of  water  or  moist  object  is  exposed  to  an  open  environment.  Under  these  conditions,  the 
water  will  begin  to  evaporate.  This  is  due  to  the  natural  tendency  of  liquid  water  to  achieve  phase 
equilibrium  with  the  moisture  content  of  the  environment  [2,  3].  As  water  evaporates,  the  latent 
heat  of  the  vaporized  water  (or  heat  of  vaporization)  is  absorbed  from  the  surrounding  environment. 
In  the  absence  of  other  mechanisms  of  heat  transfer,  a  net  cooling  effect  of  the  object's  surface  is 
experienced. 

It  is  believed  that  high  moisture  bearing  fibers  (up  to  30  wt  %  moisture)  of  polyamidobenzimidazole 
were  used  in  intercontinental  ballistic  missile  rocket  motor  casings  to  limit  thermal-induced  damage 
from  laser  threat  [4].  Morgan  et  al.  showed  that  wet  composites  can  dissipate  up  to  an  order 
magnitude  greater  thermal  energy  than  dry  composites  [5]. 

In  the  present  work,  the  feasibility  of  evaporative  cooling  of  the  surface  of  a  composite  panel  is 
numerically  assessed  as  a  means  of  reducing  the  external  temperature  of  military  aircraft.  In  the 
present  formulation,  we  develop  a  simple  one-dimensional  model  of  the  evaporative  cooling  process. 
We  consider  active  evaporative  cooling,  where  active  implies  that  the  moisture  is  actively  fed  from 
an  internal  reservoir  to  the  outer  surface  of  the  panel  by  means  of  an  active  mechanical  process. 
The  governing  differential  equation  describing  the  transient  heat  transfer  is  derived  from  the 
principle  of  energy  conservation.  A  Galerkin  weak  form  finite  element  model  is  used  to  solve  the 
resulting  nonlinear  equation.  We  present  transient  surface  temperature  results  for  a  variety  of 
environmental  conditions. 


Physical  Model 

An  illustration  of  the  physical  model  is  provided  in  Figure  1.  To  simplify  the  resulting  governing 
equation  we  consider  an  isotropic  and  homogeneous  flat  panel  or  plate  of  dimensions  L  '  L  '  L 

where  L  ?  L .  Since  the  thickness  is  much  smaller  than  the  other  dimensions  of  the  plate,  we 
will  assume  that  the  current  state  of  the  panel  is  a  function  of  the  thickness  coordinate  x  and  time  t 
only.  We  therefore  consider  an  open  one-dimensional  region  W=  (0 ,L)  with  closure  W=  | ),L|| 
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Figure  1.  The  physical  model. 


The  plate  possesses  a  uniform  initial  temperature,  T  (x,0)  =  T0.  At  t  =  0 ,  the  plate  is  exposed 

on  both  surfaces  to  streams  of  air  containing  known  free  stream  temperatures  and  moisture  contents. 
At  x  =  L ,  the  plate  is  exposed  to  an  impinging  stream  of  air  at  an  elevated  temperature  where  heat 
transfer  is  induced  by  forced  convection  of  the  impinging  fluid.  At  x  =  0,  the  panel  is  simply  exposed 
to  stagnant  air  where  convective  heat  transfer  occurs  by  buoyancy  driven  fluid  flow.  In  addition  to 
these  effects,  liquid  water  is  actively  driven  at  a  known  rate  through  the  panel  where  it  evaporates  at 
x  =  L  .  For  the  purpose  of  our  numerical  analysis,  we  assume  that  the  rate  of  moisture  flux  can  be 
mechanically  controlled.  We  further  assume  that  the  water  leaves  the  surface  (at  x  =  L  )  by 
evaporation  only.  A  more  detailed  mathematical  description  of  the  heat  transfer  process  is  provided  in 
subsequent  sections. 


Derivation  of  Governing  Heat  Equation 

In  this  section  we  present  a  derivation  of  the  governing  equation  for  the  coupled  heat  and  mass  transfer 
arising  from  the  transport  of  liquid  water  through  a  composite  panel.  We  proceed  by  developing  the 

equations  for  a  general  continuum  in  j  3  and  then  specialize  the  formulation  to  the  one-dimensional 
case.  To  this  end  we  consider  an  open  bounded  material  region  WI  j  3  with  boundary  G=  ^fW  and 
closure  W=  WUG.  The  boundary  can  be  expressed  as  G=  GD  U  G^  where  GD  I  G^  =  JE,  and 

Gfl  and  G^  are  the  Dirichlet  and  Neumann  boundary  conditions  respectively.  Applying  the 

conservation  of  mass  and  energy  principles  to  the  material  points  in  W  results  in  the  following 
expressions 


—  +  N  x(rv)=  0 
V  w 
De  ~  '■ 

r  —  =  s  :  D  -  N  xq 
Dt 


in  W  (0,  7  u 


(1) 
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In  Eq.  (l)  r  is  the  density,  v  is  the  velocity  vector,  e  is  the  internal  energy,  s  is  the  Cauchy 
stress  tensor,  D  is  the  symmetric  part  of  the  velocity  gradient  tensor  and  q  is  the  heat  flux  vector.  In 
addition,  the  time  derivative  operator  in  Eq.  (l)  is  the  material  time  derivative  operator  and  the 
divergence  operator  (i.e.,  N  x)  acts  with  respect  to  the  spatial  coordinates  x.  It  is  convenient  to  express 

the  internal  energy  in  terms  of  the  enthalpy  h  and  thermodynamic  pressure  p  as 

e  =  h  -  —  (2) 

r 


In  addition  it  is  useful  to  decompose  the  Cauchy  stress  tensor  into  the  sum  of  its  hydrostatic  and 
deviatoric  parts 

S  =  -  pi  +  t  (3) 


where  t  is  the  deviatoric  stress  tensor, 
expressed  as  £6^] 


Dh 

r - 

Dt 


Using  Eqs.  (l)i,  (2)  and  (3)  allows  the  energy  equation  to  be 
Dp 

—  =  t  :  D  -  N  xq  (4) 

Dt 


For  the  present  analysis  the  effects  of  p  and  t  will  be  ignored.  This  leads  to  the  following  simplified 
energy  equation 


r  —  +  r  v  xN/7 =  -  N  xq 
1' 


(5) 


Applying  the  principle  of  mass  conservation  to  the  liquid  water  driven  through  the  flat  panel 
results  in  the  following  standard  continuity  equation 

Ifr 

— 551  +  N  x(r  v  )  =  0  (6) 

v  w  w7  V  / 


In  Eq.  (6),  r  is  the  partial  density  of  the  liquid  water  and  Vw  is  the  velocity  of  the  same.  In  the 
present  formulation  Vw  is  taken  as  a  known  quantity.  As  a  result  it  is  unnecessary  to  solve  Eq.  (6). 


Rule  of  Mixtures  and  Constitutive  Relationships 

The  material  in  this  analysis  is  treated  as  a  binary  mixture  of  liquid  water  imbedded  in  a  solid  epoxy 
composite  plate.  Quantities  associated  with  the  bulk  material  such  as  enthalpy,  thermal  conductivity, 
velocity,  etc.,  are  expressed  throughout  this  work  using  the  rule  of  mixtures.  An  arbitrary  quantity  / 
associated  with  the  bulk  material  is  therefore  approximated  as 

rj  =  r  j  +  r  j  (7) 

where  the  subscripts  e  and  w  refer  to  the  epoxy  matrix  and  the  liquid  water,  respectively.  The  quantity 
r  in  Eq.  (7)  is  the  partial  density  of  a  particular  species.  The  partial  density  is  defined  as  the  mass  of 

constituent  i  in  a  small  representative  volume  of  the  mixture  divided  by  the  total  representative  volume 
£7].  Eq.  (7)  can  be  used  to  obtain  the  following  expression 

ryxN/i  =  r  v  xN/7  (8) 

The  constitutive  equations  employed  are  Fourier's  law  of  heat  conduction  and  the  standard 
enthalpy-temperature  relationship.  For  the  isotropic  case  these  equations  assume  the  following  forms 
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q  =  -  kN  T ,  h  =  cT  (9) 

where  k  is  the  thermal  conductivity  and  Cp  is  the  effective  specific  heat  of  the  bulk  material.  Using 
Eqs.  (8)  and  (9)  in  Eq.  (5)  results  in  the  following  energy  equation  for  the  mixture 

rc  rev  xNJ  =  N>4nt)  in  W  (O ,T$  (10) 

p  ^  w  P,w  w  V  /  V  U  v  ' 

where  Cpw  is  the  specific  heat  of  water.  In  Eq.  (10)  we  have  assumed  that  Cp  is  constant.  For  the 
present  analysis  we  will  be  concerned  with  the  one  dimensional  form  of  Eq.  (10)  which  can  be  expressed 


rc  ^-  +  r  c  v  51  =  XjllZlJ  in  (o,L)'  (0 


Finite  Element  Model  of  Energy  Equation 
Weak  Form  of  Energy  Equation 

In  this  section  we  develop  a  weak  form  Galerkin  finite  element  model  of  the  energy  equation  given  by 
Eq.  (ll).  In  our  finite  element  formulation  we  discretize  the  computational  domain  W=  |),Ljsj  into  a 

set  of  NE  non-overlapping  subdomains  (or  elements)  W  =  such  that  W=  "[J^W  .  The 

variational  or  weak  form  problem  is  to  find  T  1  H 1  (W )  for  all  test  functions  w  I  H I  (W  )  such  that 

B  (w,T  )  =  /  (w)  (12) 

where  the  bilinear  and  linear  forms  are  given  as 


B  (w,T  )=  5*‘!-c  w^-+  v  )T  \lx 

y  }  °*.  i p  v  n*  n*  w  jw w)  s  (is) 

l(W)=  W(Xa)Ql  +  W(xh)Q2 

In  the  above  expressions  we  have  made  use  of  the  standard  definition  of  the  Sobolev  space 
H 1  (W )  =  IT  12  (W )  .  The  space  H\  (w)  associated  with  the  test  functions  w  is  defined  as 

Hi  (w)=  {w  I  H 1  (W ) :  w  =  0  on  G).  where  Cf  is  the  element  boundary  where  the  temperature 
is  specified. 

For  the  sake  of  brevity  we  have  dropped  superscripts  e  from  quantities  in  the  above  equations  and 
throughout  the  remainder  of  this  work.  It  should  therefore  be  remembered  that  Eq.  (12)  applies  to  the 

eth  element  of  the  finite  element  model.  Assuming;  r  ,  c  and  V  to  be  constant  allows  the  bilinear 

7  p,w  W 

form  to  be  reduced  to 
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(14) 


X.CC 

B  ( w,T  )=  o  §rc, 


32  F  ,  ,1w1T  1™^ 

w  - —  +  k - rev  - — T  irlx 

P  m+  *Tv  «Tv  w  p,w  w 


V  lx  1 X 


lx  0 


The  quantities  Q  are  the  generalized  element  forces  which  can  be  expressed  as 


8  1 T  V 

Q.  —  -  *  - - r  c  Tv  V 

1  g  w  p,w  w  p 

e  u=x 


e  f  r  v 

Q=*k  - - r  c  Tv  V 

2  e  w  p’w  M'u 


u= 


(15) 


Semi-Discrete  Finite  Element  Equations 

In  the  present  formulation  we  assume  a  space-time  decoupled  finite  element  solution.  The  temperature 
field  in  the  weak  form  is  therefore  approximated  within  each  finite  element  using  the  following 
interpolation  scheme 


T(x,t)=  a  TJ  (x) 


(16) 


7=1 


where  xl  In  Eq.  (16),  y .  (x)  are  the  one-dimensional  Lagrange  interpolation  functions 


which  can  be  expressed  as 


'  (v)  0" 


x  -  x, 


X  -  X, 


(17) 


where  x  —  2 (x  -  x  )/  h  -  1  and  h  —  x,  -  x  .  The  quantities  x,  are  the  locations  of  the  nodes  in 

\  a  /  x  x  b  a  1  k 

terms  of  the  natural  coordinate  x  We  employ  the  standard  Galerkin  approach  and  define  a  set 

of  weighting  functions  w.  (x)  as 

w(x)=y,.(x)  (is) 

The  semidiscrete  finite  element  equations  are  obtained  by  inserting  Eqs.  (16)  and  (18)  into  Eq.  (14) 
which  results  in  the  following  element  equations 

rsP)+  ¥it}=  k'l  <i9> 

The  finite  element  matrices  appearing  in  Eq.  (19)  can  be  determined  from  the  following  formulas 

Ce  =  ’re  y  .y  dx 

ij  U  p-  j 


as 


Kl  =  6, 


dy  dy 


dy . 


dx  dx 

F,'  =  y,(*.)e,  +  y,(x>)e 


r  c  v  — Adx 


w  p,w  w  dx  J  ^ 


(20) 


The  standard  finite  element  assembly  operator  |J3]  is  used  to  construct  the  global  set  of  finite  element 
equations  which  assume  the  following  general  form 

rsM+p^}=ki  <2>) 
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Fully  Discretized  Finite  Element  Equations 


The  fully  discretized  finite  element  equations  are  obtained  by  partitioning  the  time  interval  |),  t  j|I  j 
of  interest  into  a  set  of  N  non-overlapping  subintervals  such  that 


(22) 


The  solution  is  obtained  by  solving  an  initial  value  problem  within  each  subregion  |sTs+1jj  where  the 

solution  is  known  at  .  In  seeking  a  solution  to  Eq.  (21)  within  |s,TJ+1j'j  we  replace  the  time 

derivative  operator  with  a  discrete  approximation  through  the  use  of  the  a  -family  of  time 
approximation 

(!-)«,  + oft..  <*»> 


where  0  £  a  £  1 .  Using  Eq.  (23)  in  Eq.  (21)  results  in  the  following  fully  discretized  finite  element 
equations 


+  aDt 


run., 


-  d'.(>-  a)vtM'(nn-  n) 


(24) 


In  the  present  analysis  we  use  the  Crank-Nicolson  scheme  where  a  =  1  /  2 .  For  linear  systems  this 
results  in  a  temporal  integration  scheme  that  is  unconditionally  stable  with  second  order  accuracy  with 

respect  to  the  time  step  f  -  t 

evolves  with  time.  Assuming  j|  to  be  constant  within  a  given  time  interval  yields 


In  Eq.  (24)  we  have  made  no  assumptions  on  how 


(it,  +  aD'.pUU 


'D'Ti,+fLn 

da>-  «xpsn-  xi) 


(25) 


Initial  and  Boundary  Conditions 


To  complete  the  finite  element  formulation  it  is  necessary  to  specify  appropriate  initial  and  boundary 
conditions.  For  the  initial  conditions,  we  assume  that  the  temperature  in  the  plate  is  given  as 


T  (t,o)=  r0  *!jug 


(26) 


where  T  is  a  real  constant. 

For  the  boundary  conditions,  we  consider  convection  driven  heat  transfer  from  both  top  (x  =  L  ) 
and  bottom  (x  =  0)  surfaces  of  the  flat  panel.  The  mass  flux  of  the  water  through  the  plate  thickness 
is  considered  uniform  at  a  given  moment  in  time.  This  flux  is  required  to  be  equal  to  the  moisture 
evaporation  flux  at  the  top  surface  of  the  plate.  The  mass  flux  of  the  moisture  is  determined  by  the 
following  formula 

=  r  v  =  bh  (r  -  r  v  )  (27) 

w  w  w  m  \  w,#  /  ^  / 

where  h  is  the  mass  transfer  coefficient  and  r  is  the  saturated  vapor  density  of  the  evaporating; 

m  w,5  i  j  i  o 


11 


fluid.  The  quantity  r  is  the  partial  density  of  water  vapor  in  the  air  surrounding  the  panel.  It  is 
important  to  note  that  this  quantity  is  a  function  of  the  surface  temperature  of  the  plate.  The 
parameter  b  =  ^),  1^  is  the  ratio  of  the  “wetted”  area  of  the  plate  surface  (at  x  =  L  )  to  the  actual 

surface  area.  When  b  =  1  the  surface  is  considered  to  be  fully  saturated.  Using  Eq.  (27)  allows  the 
element  stiffness  matrices  to  be  expressed  as 


Ke 

ij 


O 


dx  dx 


-  bh 


w,¥ 


\  dy  | 

)c  — -y  mix 

/  p,w  I  ^  — 


dx  '  1  -k 


(28) 


The  boundary  conditions  for  the  finite  element  equations  can  be  expressed  in  the  following  forms 


e;  =  i(r,  -  T)+  bhm{ 


r  -  r 

w,.?  w,  ¥ 


)c  rj 

/  />.»  g 


(29) 


Q 


NE  _ 
1 


w,¥ 


+  c  T 

p,w 


(30) 


where  h  is  the  convection  heat  transfer  coefficient  and  T  is  the  farfield  temperature  of  the  air.  It  is 

important  to  note  that  the  values  of  h  and  T  appearing  in  Eq.  (29)  are  not  necessarily  the  same  as 

those  appearing  in  Eq.  (30)  (see  Figure  1).  The  boundary  conditions  are  incorporated  into  the 
assembled  finite  element  equations  by  applying  the  following  modifications  to  the  global  stiffness  matrix 
and  force  vector 


K  =  K  +  h  -  bh  ( r  -  r  )c 

11  11  m  \  w,S  w,¥  /  p,W 

K  =  K  +  h  +  bh  (r  -  r  )c 

NN  NN  m  \  w,s  vf,¥  /  p, i 


Fx  =  hT ¥ 

FNN  -  hTZ  -  bK  (rw,-  -  rw,¥  )hfg  (32) 

where  NN  is  the  total  number  of  nodes  in  the  finite  element  model. 

Numerical  Results  and  Discussion 

In  this  section  we  present  numerical  results  obtained  from  finite  element  simulations  conducted  using 
the  model  described  in  the  previous  section.  The  finite  element  model  is  implemented  using  the 
MATLAB  programming  language.  Double  precession  accuracy  is  used  for  all  quantities  in  the  model. 
For  all  simulations  conducted  we  have  employed  a  mesh  containing  25  equally  spaced  cubic  elements 
(p  =  3).  This  amounts  to  101  total  degrees  of  freedom  (or  nodes)  in  the  model.  This  fine  mesh  size 
is  necessary  in  the  weak  form  Galerkin  formulation  to  avoid  spurious  oscillations  in  the  numerical 
solution,  due  to  the  advection  term  that  is  present  in  the  energy  equation.  The  solution  procedure  is 
iterative  since  the  saturation  density  of  water  vapor  is  a  nonlinear  function  of  the  surface  temperature. 
Other  nonlinearities  are  also  present  in  the  model,  as  other  material  parameters  are  also  functions  of 
temperature.  Curve  fits  of  these  parameters  are  presented  in  the  Appendix  of  this  report. 

Table  1.  Dimensions,  parameters  and  material  properties  for  finite  element  model. 

Parameter  Description 

L  =  0.01  Panel  thickness  [)nT] 

L  —  0.06  Length  of  plate  [nr] 

p 

H  =  0.051  Distance  from  exit  of  nozzle  to  flat  plate  ]m] 
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d,  =  0.009525 

h 

r  =  0.06 
c  =  1,000 

p,e 

c  -4,181 

p,w 

r  =  1,130 

e 

k  =  0.16 

e 

V  =  9.935 
g  =  9.81 


Inner  diameter  of  nozzle  jjnj 

Radial  distance  for  averaging  of  h  [m] 

Specific  heat  of  epoxy  £J/(kg  K)[] 

Specific  heat  of  liquid  water  [21/ (kg  K)] 

Density  of  epoxy  panel  £kg/ m3^] 

Tliermal  conductivity  of  epoxy  [WV(m  K)(] 

Freestream  velocity  of  air  impinging  against  panel  \jn/ sj] 
Gravitational  constant  [hn/ s2[] 


In  addition  to  the  parameters  given  in  Table  1,  we  specify  the  freestream  air  temperature  of  the 
impinging  jet  to  be  given  as  T  =  90° C  .  For  the  numerical  results,  we  perform  a  parametric  study  of 
the  evolution  of  the  surface  temperature  of  the  panel  by  varying  b  and  the  relative  humidity  /  of 
the  freestream  impinging  air.  The  relative  humidity  /  is  the  ratio  of  the  mass  of  the  moisture  in  the 
air  to  the  maximum  amount  of  moisture  the  air  could  hold  at  the  same  temperature. 


<i>  =  1.00 
<|>  =  0.75 
<|>  =  0.50 
<|>  =  0.25 
<|>  =  0.00 


Figure  2.  Surface  temperature  profiles  (at  x  =  L)  of  an  epoxy  panel  exposed  to  an  impinging  jet  [T  -  25°C, 
r¥  =  90° C  and  b  =  1.00). 


In  Figures  through  5  we  present  transient  surface  temperature  results  for  flat  panels  that  each 
posses  an  initial  temperature  TQ  =  25°C.  A  constant  time  step  of  Dr  =  0.05  s  is  employed  in  all 

cases.  In  each  figure  a  different  surface  saturation  ratio  (i.e.,  b  =  |). 25, 0.50, 0.75, 1.00^  is  evaluated. 

Furthermore  in  each  figure  we  evaluate  the  effect  of  the  relative  humidity  /  of  the  environment  on 
the  overall  surface  temperatures.  As  expected,  the  effects  of  evaporation  play  a  major  role  in  reducing 
the  surface  temperature  profiles  for  all  cases.  It  can  be  seen  that  under  ideal  circumstances,  where  the 
moisture  at  the  surface  of  the  panel  is  maintained  at  saturation  conditions  (i.e.,  cases  where  b  =  1.00), 
significant  cooling  of  the  surface  is  experienced.  When  the  freestream  air  is  dry,  evaporation  maintains 
the  surface  temperature  very  close  to  the  initial  temperature  of  the  plate. 
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<]>  =  1.00 
<|>  =  0.75 
<]>  =  0.50 
<]>  =  0.25 
<|>  =  0.00 


Figure  3.  Surface  temperature  profiles  (at  x  =  L)  of  an  epoxy  panel  exposed  to  an  impinging  jet  (T  —  25°C, 
r¥  -  90° C  and  b  =  0.75 ). 


For  less  ideal  cases  (i.e.,  cases  where  b  <  1.00)  the  cooling  effects  are  less  pronounced,  but  still 
significant  as  can  be  seen  in  Figures  ,  and  .  We  observe  a  similar  trend  when  considering  the 
influence  of  the  relative  humidity  of  the  external  flow.  More  humid  environments  lead  to  less 
evaporation  from  the  panel  surface,  and  ultimately  lower  cooling  effects. 


<j>  =  1.00 
<|>  =  0.75 
<|>  =  0.50 
<|>  =  0.25 
<|>  =  0.00 


Figure  4.  Surface  temperature  profiles  (at  x  =  L)  of  an  epoxy  panel  exposed  to  an  impinging  jet  (T  —  25°C , 
r¥  =  90° C  and  b  =  0.50). 
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<]>  =  1.00 
<|>  =  0.75 
<]>  =  0.50 
<]>  =  0.25 
<|>  =  0.00 


Figure  5.  Surface  temperature  profiles  (at  x  —  L)  of  an  epoxy  panel  exposed  to  an  impinging  jet  (T  —  25°C , 
r¥  -  90° C  and  b  =  0.25 ). 

It  is  also  important  to  quantify  the  mass  of  the  moisture  lost  to  the  environment  as  a  result  of 
the  evaporative  cooling  process.  In  Table  2  below  we  present  values  for  the  total  mass  (per  unit  area) 
of  water  lost  due  to  evaporation  during  the  numerical  simulations  as  a  function  of  the  panel  wetness  b  , 
and  relative  humidity  /  of  the  environment.  Approximate  average  mass  flux  rates  can  be  obtained  by 
dividing  the  values  in  Table  2  by  the  total  simulation  time  (t  =  500  s). 

Table  2.  Effect  of  panel  wetness  b  ,  and  relative  environmental  humidity  /  ,  on  total  mass  (per  unit  area)  of 
water  vaporized  in  evaporation  process  after  500  s. 

Mass  of  water  evaporated  after  500  s  (kg/m2) 


f  b  =  1.00  b  =  0.75  b  =  0.50  b  =  0.25 


1.000 

0.000 

0.000 

0.000 

0.000 

0.750 

0.490 

0.434 

0.359 

0.246 

0.500 

0.631 

0.571 

0.490 

0.359 

0.250 

0.716 

0.655 

0.571 

0.434 

0.000 

0.776 

0.716 

0.631 

0.490 

In  appraising  the  numerical  results  presented  in  this  study  it  is  important  to  draw  attention  to 
two  of  the  key  assumptions  inherent  in  the  numerical  model.  First;  in  the  present  study,  the  convective 
mass  transfer  coefficient  is  determined  using  the  well-known  analogy  between  convective  heat  and  mass 
transfer  as  discussed  in  the  Appendix  of  this  report.  Actual  mass  transfer  coefficients  may  deviate  to 
some  degree  from  the  values  obtained  using  the  heat  and  mass  transfer  analogy.  The  second 
assumption  is  that  the  plate  surface  remains  saturated  and  the  water  leaves  the  surface  only  by 
evaporation.  In  other  words,  we  have  assumed  that  the  water  is  not  blown  from  the  surface  by  the 
impinging  jet  prior  to  evaporation. 

Conclusions 

In  this  work  we  have  presented  a  numerical  model  to  assess  the  potential  of  active  evaporative  cooling  as 
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means  of  reducing  the  external  temperature  of  military  aircraft.  We  introduced  a  simplified  physical 
model  of  the  phenomena  consisting  of  a  flat  panel  exposed  to  an  impinging  jet  of  hot  air  on  one  surface. 
Evaporative  cooling  was  induced  in  the  theoretical  model  by  actively  driving  moisture  through  the  panel 
to  the  outer  surface.  The  governing  heat  equation  for  the  physical  model  was  derived  from  the 
conservation  of  energy  principle  and  appropriate  constitutive  relationships.  The  equation  was  solved 
numerically  using  a  weak  form  finite  element  formulation  in  conjunction  with  the  a  -family  of  time 
approximation  and  the  Crank-Nicolson  time  integration  scheme. 

Parametric  studies  were  conducted  on  a  benchmark  problem,  where  the  panel  wetness  b  ,  and 
relative  environmental  humidity  /  were  varied.  Surface  temperature  profiles  were  then  plotted  to 
demonstrate  sensitivity  of  the  panel’s  surface  temperature  to  these  parameters.  The  results  indicate 
that  under  ideal  situations,  where  there  was  low  relative  environmental  humidity,  and  panel  surface 
saturation  near  unity,  active  evaporative  cooling  is  able  to  greatly  reduce  a  panel's  surface  temperature. 
However,  we  can  also  conclude  that  under  less  ideal  circumstances,  the  evaporation  process  still  counters 
much  of  the  heat  of  the  impinging  fluid. 
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MICROCHANNEL  MODEL  EXPERIMENTS 

Many  problems  with  thermal  management,  signature  reduction,  and  self-repairing  materials  require 
micro-channels  in  the  structure.  These  micro-channels  would  allow  fluids  to  flow  through  the 
structure  to  provide  repair  agents,  collect  waste  heat,  and  to  reduce  surface  temperature.  Our 
experiments  attempted  to  create  a  model  system  of  microchannels. 

It  is  hoped  that  various  directions  of  fluid  flow  through  complex  microchannel  networks  will 
conduct  heat  reduce  temperatures.  Such  system  might  dissipate  enough  energywith  effective  heat 
transfer  through  a  large  surface  area  to  volume  ratio.  Furthermore,  high  thermal  conductivity  of 
carbon  fiber  might  enhance  heat  conduction  to  the  microchannel  network. 

Sample  preparation 

MicroChannel  experiments  were  conducted  with  these  materials: 

•  Sample  material  (epoxy/T-300  carbon  fiber) 

-  Sample  size  :  76.2  mm  *  76.2  mm  *  2  mm  (L*W*H) 

-  Diameter  of  micro-channel :  100  ~  220  pm 

-  Sample  fiber  direction  :  90° 

•  Sample  preparation  process  (micro-channel  fabrication) 

-  Putting  lubricated  nickel  chromium  wires  through  the  sheets  of  composite  positioning  fiber 
direction  before  curing 

-  Applying  electric  current  both  end  of  wire  to  burn  the  material  around  the  wire  which 
makes  the  wire  pull  out  smoothly. 

-  Pulling  the  wires  from  the  composite  to  get  micro-channel 
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Figure  6  (Fig  1).  Nickel  chromium  wires  in  a  layered  composite 


Figure  7.  Variable  electric  current  generation 


Figure  8.  Epoxy/T-300  carbon  fiber  with  micro-channel 

Figure  6.  shows  a  prepared  sample  composite  material  with  wires  after  curing.  To  release  the  the 
wires  smoothly,  electric  current  was  applied  through  the  wires.  Figure  7.  shows  electric  current 
control  with  a  variac.  Finally,  microchannels  in  the  composite  were  fabricated  as  Figure  8  shows. 
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Complex  microchannel 


First,  90°  fiber  angles  between  layers  provides  90°  microchannel  direction  changes  from  layer  to 
layer.  Microchannels  positioned  on  the  adjacent  sheets  might  meet  and  mix  fluid  internally. 


Heat  Source 


Figure  9  .  Schematic  diagram  of  90°  internal  micro-channels 

Mathematical  modeling  of  conductive  heat  transfer  in  composites 

In  this  section,  a  mathematical  modeling  of  heat  transfer  in  composite  is  introduced.  Composite 
laminate  multiple  fiber  angles  and  analyzing  the  thermal  properties  of  oriented  fibrous  composite  is 
difficult,  since  fiber  orientation  affects  thermal  properties.  Each  lamina  is  transversely  isotropic 
about  its  fiber  axis.  Consider  that  this  case  is  as  two-dimensional  heat  transfer,  then  /e,yis  the  thermal 
conductivity  of  any  lamina  [1]. 

kij  =  (ka  +  kt)  piPj  +  ktSij  (1) 

ka :  axial  heat  conductivities 
kt\  transverse  heat  conductivities 
Pi :  a  unit  vector  to  axial  fiber  direction 
5ij\  Kronecker  delta 

From  the  mathematical  model,  kt  (transverse  heat  conductivities)  is  an  important  property  since 
heat  source  from  the  surface  area  transfer  to  the  media  transversly  [2]. 

To  find  n-th  layer  heat  conductivities, 

kitx=  X  knxxhn/  £  hn  (2) 
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kxx :  average  of  x  direction  heat  conductivities 
knxx :  the  n-th  layer  heat  conductivity 
hn :  thickness 

Fluid  flow  may  be  very  slow  and  should  be  a  laminar  flow  since  cross  sectional  area  of  the  channels 
is  micro  size.  As  a  result  of  limitation  of  fluid  velocity  conductive  heat  transfer  dominates  in  the 
channels.  Therefore,  capillary  action  is  expected,  but  microchannels  in  each  layer  meet  neighboring 
microchannels  in  90°  and  are  expected  to  mix  fluid  actively  which  makes  convective  heat  transfer 
increase.  Furthermore,  if  fiber  angles  more  complicated  than  90°,  more  enhancement  of  heat 
transfer  will  be  expected. 
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APPENDIX 


A.  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS 

The  convective  heat  transfer  coefficients  li  are  determined  from  a  combination  of  analytical  and 
empirical  formulas  that  are  readily  available  in  the  literature  [2, 7,9].  In  the  case  of  a  single  round 
nozzle  discharging  fluid  against  a  flat  surface  (i.e.,  a  single  impinging  jet),  the  Reynolds  and  average 
Nusselt  numbers  are  defined  as  []7[] 

Vd. 

Re  °  — (A.l) 


—  hd, 
Nu  °  — ^ 


(A.2) 


where  dh  is  the  inner  diameter  of  the  nozzle  from  which  the  gas  stream  is  flowing,  k  is  the  thermal 
conductivity  of  the  air  and  V  is  the  air  velocity  The  following  expression  is  obtained  from  the 
literature  for  approximating  the  average  convection  heat  transfer  coefficient  h  []7[] 

2ka  Re1/2  Pr0'42  (l  -  l.k/;i/  r)(l  +  0.005  Re0'55  ^ 


,1/2 


h  = 


+ 


0-1  6 >/„/,'* 


(A.s) 


where  the  Prandtl  number  is  defined  in  general  as 

v  r  c 

pr  o  a  a  p 


(A.4) 


In  Eq.  (A.4),  the  quantities  v j and  c  are  the  kinematic  viscosity,  density  and  specific  heat  of  air. 

Eq.  (A.S)  is  used  to  estimate  the  average  convection  heat  transfer  coefficient  h  in  the  vicinity  of  the 
stagnation  point  of  the  flow.  The  quantity  r  is  the  radial  distance  from  the  stagnation  point  to  an 
arbitrary  distance  on  the  surface  of  the  plate  over  which  h  is  averaged.  The  quantity  H  is  the 
distance  from  the  exit  plane  of  the  nozzle  to  the  surface  of  the  flat  plate. 

It  is  also  necessary  to  determine  the  average  convective  heat  transfer  coefficient  for  the  bottom 
side  of  the  plate  (x  =  0)  where  convection  occurs  due  to  buoyancy  driven  flow  conditions.  For  free 
convection  (i.e.,  buoyancy  driven  flows)  the  averaged  Nusselt  number  is  given  as  [To)] 

—  hL 

Nu  °  (A .5) 

k 


In  addition,  the  Rayleigh  number,  Rcil  ,  is  defined  as  [To]] 


Rcil 


gr 


c  \r  - 

a  p,a  s 


L 


vkT 

a  a  j 


(A.6) 


where  gis  the  gravitational  constant.  The  averaged  Nusselt  number  is  given  as  [To]] 


Nu  =  0.68  + 


0.6707?a 


1/4 


9/16u/9 


RaL  <  109 


(A.7) 


#+  (0.492/  Pr)  u 

g  v  ;  a 
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The  average  convective  heat  transfer  coefficient  h  for  free  convection  may  therefore  be  determined 
from  Eqs.  (A.4)  through  (A.7). 


B.  CONVECTIVE  MASS  TRANSFER  COEFFICIENTS 

It  is  also  necessary  to  determine  the  convective  mass  transfer  coefficient  h  .  To  do  so  it  is  important 
to  introduce  some  non-dimensional  parameters.  The  Schmidt  number  Sc  is  defined  as  £7^] 

Sc  °  (B.l) 

D 


For  parallel  flow  over  a  flat  plate  the  average  Sherwood  number  is  given  as  £7] 

—  h  L 

Sh°  m  p 


D 


(B.2) 


where  the  quantity  D  is  the  diffusivity  of  water  vapor  in  air.  In  the  case  of  a  single  impinging  jet, 
the  average  Sherwood  number  is  expressed  as  ^7^j 


—  h  d, 
Sh  °  m  h 


D 


(B.3) 


There  are  strong  analogies  between  heat  and  mass  transfer.  This  often  leads  to  a  similarity 
between  the  convective  heat  and  mass  transfer  coefficients.  This  correspondence  is  of  the  form  £2,7,9)] 

Nu  Sh 

(B.4) 


/V 


Scn 


where  Nu  and  Sh  may  be  local  or  averaged  quantities.  If  the  heat  transfer  coefficient  is  known,  the 
mass  transfer  coefficient  may  be  determined  from  Eq.  (B.4).  The  above  analogy  holds  for  low  rates  of 
mass  transfer  £2, 9]].  Eq.  (B.4)  can  be  used  to  obtain  the  following  expression  for  the  convective  mass 
transfer  coefficients 


h  = 


h 


..1-  } 

3d  k  2. 

r  c  £ - T 

fl  p'a  gr  c  D  * 

®  a  p,a  a  Vj 


(B.5) 


For  parallel  flow  over  a  flat  plate  it  has  been  shown  that  n  —  1/  3  [)9^j.  For  the  case  of  flow  from 
impinging  jets  it  has  been  suggested  that  n  =  0.42  [)9^j. 


C.  MATERIAL  PROPERTIES 

Many  of  the  material  properties  and  coefficients  used  in  this  work  are  not  constant.  Some  vary 
significantly  with  temperature  and/ or  pressure.  In  this  work  variables  are  assumed  to  be  a  function  of 
temperature  only.  The  following  are  curve  fits  of  some  of  the  published  properties  of  air  [)7j]  at 
atmospheric  pressure 
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r  =-  7.9282'  10'  147; 5  +  9.5330'  10  '  n7  4  -  5.0164'  10'8'73 

a  f  f  f 

+  1.68472'  10'  57  2  -  0.00477  +  1.2770 
c  =-  8.8889'  10147  6  +  9.9528'  10’  nT; 5  -  4.2892'  10'  87 4 

f  f  f 

+  8.2319'  10‘ *7/-  0.00027^+  0.01537^+  1006.5729 
v  =-  4.1026'  10'  197/  +  3.1179'  10' 167,4  -  1.3269'  10' 137f3 

a  f  f  f 

+  1.2036'  10'107/2  +  8.8727'  10'^  +  1.3414'  10' 5 
k  =  1.8461'  10"  157.5  -  1.4637'  10' 127,4  +  3.6606 '  10' 107,3 

a  f  f  f 

-  5.7872'  10' 187f2  +  7.93825'  10'  57;  +  0.0242 


(C.l) 


(C.2) 


(C.3) 


(C.4) 


Pr=  1.3778'  10'  'sTf6  -  1.2822'  10' 127f5  +  4.3350 '  10'  107  4 

-  6.31081'  10' 8773  +  3.8276'  1  O'  (,Tf2  -  2.3383'  10' 47^  +  0.71 17 


(C.5) 


Eqs.  (C.l)  through  (C.5)  are  valid  between  -23C  and  300C.  The  film  temperature  T f  —  —(T  +  7¥  ) 
possesses  units  of  Celsius.  Tlie  diffusion  coefficient  of  water  vapor  in  air  is  taken  from  the  literature 

DO  as 

.2.072 


D  =  1.87'  10" 10  (7  +  273)" 


(C.6) 


The  latent  heat  of  vaporization  of  water  h  and  the  saturation  density  of  water  vapor  are  given 

by  the  following  curve  fits  of  the  published  data  DO 

h  =  2.750'  10'  107,6  -  2.9275'  10'  7T/  +  7.67128  '  10' 57  4 

Jg  J  J  J  /p  ?\ 

-  0.02117^+  0.900097^-  22,372.60577^  +  2,501,365.4834  1  '  ’ 

r  =  1.6978'  10'  147,.6  -  1.8031'  10' 127f5  +  4.3235 '  10'9'  74 

"V  f  f  f 

-  4.5824'  10'8'  7/3  +  1.8577  '  10'  5Tf2  +  2.1331'  lO'D,  (C.8) 

+  5.1586'  10' 3 
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